

function f = pdf_psi(wb, z, h, WB, Z)

% Returns the value of the conditional probability density function evaluated at wb given Z = z 
% Nonparametric estimation with a smoothing parameter h


%******************* Input Variables **************************************

% wb : value of the variable of interest, where the c.d.f. will be
% evaluated
% z : value of the conditional variable
% h : smoothing parameter
% WB : vector storing the outcomes of the variable of interest
% Z : vector storing the outcomes of the covariate

%******************* Output Variable **************************************

% f : value of the conditional probability density function

L = length(WB);
pz = size(Z,2);

% ************** kernels and auxiliary sums ******************************

numer = 0; % Auxiliary for cummulative sum
denom = 0; % Auxiliary for cummulative sum
% K = @(u) (1/sqrt(2*pi))*exp(-0.5*(u.^2));  % Gaussian
% K = @(u) ((35/32)^2)*((1-u^2)^3)*(abs(u)<=1); % Triweight
K = @(u) (3/4)*(1 - u.^2)*(abs(u)<=1); % Epanechnikov

for ell = 1:L
    for k=1:pz
       KF(k) = K((z(k) - Z(ell,k))/h(k+1)); % Computation of kernel function for each variable in Z
    end
       KP = prod(KF,2); % Product of kernel functions
       numer = numer + K((wb - WB(ell))/h(1))*KP; % cummulative sum for numerator
       denom = denom + KP; % cummulative sum for denominator 
       
end

hp = prod(h); % product of variables associated bandwidths
hz = prod(h(:,2:pz+1),2); % product of variables associated bandwidths (excluding dependent variable)
num=(1/(L*hp))*numer; % numerator
den=(1/(L*hz))*denom; % denominator

f = num/den; % value of density function 

end